Embedded Laplace: find the mode when the first Newton step overshoots, and throw for a non-finite log marginal - #3473

Fixes #3469. The diff is only optically large: most of it is a block of code moved into a helper and re-indentation, so hiding whitespace in the diff view helps. Stacked on #3465 (its commit shows in the diff until it merges); suggested merge order: #3438, #3465, this PR.

AI use disclosure: This was developed with the help of Claude Code (Opus 5.5, Fable 5.1).

Summary

laplace_marginal_density_est could return the Laplace approximation evaluated at the initial guess as the log marginal, or +inf / NaN, without a message (#3469). Three changes:

1. wolfe_line_search no longer accepts the zero step. The fallback best starts at step 0, which satisfies Armijo trivially, so a failed search was reported as an accepted step and the Newton loop "converged" at the starting point. The fallback now requires best.alpha() > 0.
2. run_newton_loop searches shorter steps when no step improves the objective. From θ = 0, W can be many orders of magnitude smaller than at the mode, so the Newton step $p$ overshoots by more than $1/\alpha_{\mathrm{min}}$ (min_alpha = 1e-8) and every trial step lowers the objective or overflows. The loop now retries the search with $p$ scaled by $\alpha_{\mathrm{min}}$. For a concave objective $f$, steps that short gain at most
   $$f(a + t p) - f(a) \le t \nabla f(a)^\top p \le \alpha_{\mathrm{min}} \nabla f(a)^\top p, \qquad 0 < t \le \alpha_{\mathrm{min}},$$
   so the retry stops when this bound is below tolerance or below the rounding error of $f$, or is not finite, or min_alpha ≥ 1.
3. A singular $I + \Sigma W$ throws std::domain_error. Mathematically $\det(I + W^{1/2} \Sigma W^{1/2}) \ge 1$, but when $\Sigma$ is numerically rank deficient and $\Sigma W$ swamps the identity, the matrix is singular in floating point: in the example at α = 1e9 the jitter is below the ulp of α², so all rows of $I + \Sigma W$ are equal. Then the Newton system has no solution and the LU factorization of solver 3 has zero pivots, so its log determinant is −∞ and develop returned +inf (NaN with a negative pivot). Solvers 1 and 2 already throw when their Cholesky fails in this situation; solver 3 now throws at the first factorization with a zero pivot, instead of iterating on a meaningless step and reporting the result at the end. A check after the solver chain also throws for any other non-finite log marginal (for example a non-finite log likelihood at the final iterate); it sits after the chain so that the chain cannot catch it and fall through into the next solver. Return the final Newton iterate as the Laplace latent mean, and fix the solver-3 log determinant #3438 makes the LU log determinant throw for a non-positive determinant at the end of the loop; the two merge cleanly and complement each other.

Example from the issue (y = (3, 1, 4, 1, 5), K = α² J + 1e-8 I); the reference values are Laplace approximations at the mode in 80-digit arithmetic (mpmath):

The develop values at μ = −30 are the approximation evaluated at θ = 0. At α = 1e9, $I + \Sigma W$ is singular in double precision at the mode (the jitter is below the ulp of α²), so throwing is the correct outcome. The remaining 0.034 error at α = 1e7 is not addressed here. On the N = 100 Poisson GP model where this was found (ρ = 3e27, μ = −51.6), α = 1e4 gives −311.555 instead of −21782.9, the value at θ = 0; from α = 1e8 on, where develop returned that value, +inf or NaN, the result now throws.

The line search of run_newton_loop now lives in internal::line_search_fails, so the retry is a short loop, and the solver chain is wrapped in a lambda so that the check runs after it. An alternative to the lambda, if preferred: split laplace_marginal_density_est into an _impl holding the chain and a thin wrapper holding the check.

Tests

Three new tests, all failing on develop:

• laplace_marginal_poisson_log_lpmf_newton.step_overshoots_mode compares μ = −30, α = 3e4 and 1e5 with the mpmath references (tolerance 1e-4, since the loop stops within its tolerance of the mode, where the log-determinant term still moves the value by up to 1e-5). Develop is off by about 408 in both. α = 1e5 needs change 2; α = 3e4 needs changes 1 and 2 together.
• laplace_marginal_poisson_log_lpmf_newton.singular_log_determinant expects std::domain_error at μ = 1, α = 1e9 with solver 3 and no fallthrough, where all rows of $I + \Sigma W$ are equal and the factorization has zero pivots. Develop returns +inf.
• WolfeLineSearch.NoImprovingStepFailsInsteadOfZeroStep scripts an update where every trial step lowers the objective and expects accept_ == false with status Fail. Develop returns Armijo with accept_ == true, the zero step.

All 18 files in test/unit/math/laplace pass at O=3; the header checks of both changed headers, cpplint and clang-format pass.

Side Effects

Well-conditioned problems are unaffected: where the line search succeeds, values and cost equal develop's. In ill-conditioned cases the number of likelihood evaluations can go up as well as down: up where the first Newton step overshoots and develop stopped at the starting point, since the retry searches on to the mode; down where $I + \Sigma W$ is singular, since solver 3 now fails at the first factorization.

Counts are likelihood calls (WSL, g++, -O0); otherwise the change adds scalar comparisons only, and memory is unchanged. Rows where develop's line search failed or the branch throws (all others are equal):

Behaviour changes:

• a failed line search no longer counts as convergence (same returned iterate, one factorization pass fewer);
• a singular $I + \Sigma W$ at any iteration and any non-finite log marginal throw instead of being returned, also in the laplace_latent_* functions (before Return the final Newton iterate as the Laplace latent mean, and fix the solver-3 log determinant #3438 also for a negative LU pivot of a non-singular B);
• with max_steps_line_search = 0, μ = −30 and α = 1e5 or 1e7 throw instead of returning the value at θ = 0: the shorter step found by the retry is accepted blindly and lands where the factorization is singular.

Release notes

Embedded Laplace now finds the mode when the first Newton step overshoots, and throws instead of returning +inf or NaN.

Checklist

• Copyright holder: Jáchym Barvínek jachymb@gmail.com

  The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
  - Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
  - Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

• the basic tests are passing

  • unit tests pass (to run, use: ./runTests.py test/unit)
  • header checks pass, (make test-headers)
  • dependencies checks pass, (make test-math-dependencies)
  • docs build, (make doxygen)
  • code passes the built in C++ standards checks (make cpplint)

• the code is written in idiomatic C++ and changes are documented in the doxygen

• the new changes are tested